#Project: TV polarization
#Author: Josh Ryan and Jeff Lyons
#Working Date: Winter 23
#Task: Use Imai, Kim, Wang PanelMatch analysis--analysis using ps match--for legislator level, Figure 9
#Related files: leg_level_forstata_v6_102524.dta is loaded
#Updated for PanelMatch 3.1.1

library(foreign)
library(effects)
library(sciplot)
library(sandwich)
library(lmtest)   
library(BMS)
library(rJava)
library(lattice)
library(plyr)
library(lme4)
library(gmodels)
library(MASS)
library(stringr)
library("XLConnect")
library(reshape2)
library(xlsx)
library(dplyr)
library(doBy)
library(Hmisc)
library(stargazer)
library(xtable)
library(dplyr)
library(xtable)
library(foreign)
library(dotwhisker)#https://cran.r-project.org/web/packages/dotwhisker/vignettes/dotwhisker-vignette.html
library(berryFunctions) #for the insert row function
library(ggplot2)
library(gridExtra)
library(wnominate)
library(tidyr)
library(wnominate)
library(oc)
library(haven)
library(PanelMatch)

##Figure 8 and Figure I5 (at bottom)
dat<-read_dta("leg_level_forstata_v6_102524.dta")
#Have to keep making as dataframe because of the way Panelmatch works
dat<-as.data.frame(dat)

#keep only relevant variables, doesn't matter but panelmatch throws an error
dat.np<-dat[c("np_score", "treatment", "majority", "party_id", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro", "year", "distid1")]

#have to constantly be making these as integers, super annoying
dat.np<-as.data.frame(dat.np)
dat.np$distid1<-as.integer(dat.np$distid1)
dat.np$year<-as.integer(dat.np$year)

dat.np0<-subset(dat.np, party_id==0)
dat.np1<-subset(dat.np, party_id==1)

######np score Dems

# Create PanelData object for Democrats
panel_data_dem <- PanelData(
  panel.data = dat.np0,
  unit.id = "distid1",
  time.id = "year",
  treatment = "treatment",
  outcome = "np_score"
)

# Updated for PanelMatch 3.1.1
results.mah.ideo.dem <- PanelMatch(
  panel.data = panel_data_dem,
  lag = 4,
  refinement.method = "mahalanobis",
  match.missing = TRUE,
  covs.formula = ~ I(lag(majority, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
    I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4)) +
    I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
    I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
  size.match = 10,
  qoi = "att",
  lead = 0:4,
  forbid.treatment.reversal = TRUE
)

mset.ideo.dem <- results.mah.ideo.dem$att
PE.results.ideo.dem <- PanelEstimate(sets = results.mah.ideo.dem, panel.data = panel_data_dem)
summary(PE.results.ideo.dem)
PE.results.ideo.dem<-summary(PE.results.ideo.dem) #save as dataframe
names(PE.results.ideo.dem)
coef<-PE.results.ideo.dem[ , 1] #get time coefficients
se<-PE.results.ideo.dem[ , 2] #get standard errors
lower<-PE.results.ideo.dem[ , 3] #lower confidence bound
upper<-PE.results.ideo.dem[ , 4] #upper confidence bound
time<-c(0,1,2,3,4)

results.ideo.dem<-as.data.frame(cbind(coef,se,lower,upper,time))
results.ideo.dem

plot.ideo.dem<-ggplot(results.ideo.dem, aes(time,coef)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) +
  theme_classic() +
  scale_y_continuous(breaks=seq(-.2, .2, .04)) +
  coord_cartesian(ylim=c(-.2, .2)) +
  scale_x_continuous(breaks=seq(0, 4, 1)) +
  theme(panel.border = element_blank(), panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"))+
  theme(axis.line.x = element_line(color="black", size = .5),
        axis.line.y = element_line(color="black", size = .5))+
  geom_hline(yintercept=0, color="gray")+
  xlab("Time From Adoption (in Years)") + ylab("Estimated Coefficient with CI")+
  ggtitle("Democratic Ideology")+theme(plot.title = element_text(hjust = 0.5))

plot(plot.ideo.dem)

######################GOP
# Create PanelData object for Republicans
panel_data_gop <- PanelData(
  panel.data = dat.np1,
  unit.id = "distid1",
  time.id = "year",
  treatment = "treatment",
  outcome = "np_score"
)

#covariate balance propensity score matching
results.cbps.ideo.gop <- PanelMatch(
  panel.data = panel_data_gop,
  lag = 4,
  refinement.method = "CBPS.match",
  match.missing = TRUE,
  covs.formula = ~ I(lag(majority, 1:4)) + I(lag(veto, 1:4)) + I(lag(mds1, 1:4)) + 
    I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4)) +
    I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
    I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
  size.match = 10,
  qoi = "att",
  lead = 0:4,
  forbid.treatment.reversal = TRUE
)

mset.ideo.gop <- results.cbps.ideo.gop$att
PE.results.ideo.gop <- PanelEstimate(sets = results.cbps.ideo.gop, panel.data = panel_data_gop)
summary(PE.results.ideo.gop)
PE.results.ideo.gop<-summary(PE.results.ideo.gop) #save as dataframe
names(PE.results.ideo.gop)
coef<-PE.results.ideo.gop[ , 1] #get time coefficients
se<-PE.results.ideo.gop[ , 2] #get standard errors
lower<-PE.results.ideo.gop[ , 3] #lower confidence bound
upper<-PE.results.ideo.gop[ , 4] #upper confidence bound
time<-c(0,1,2,3,4)

results.ideo.gop<-as.data.frame(cbind(coef,se,lower,upper,time))
results.ideo.gop

plot.ideo.gop<-ggplot(results.ideo.gop, aes(time,coef)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) +
  theme_classic() +
  scale_y_continuous(breaks=seq(-.2, .2, .04)) +
  coord_cartesian(ylim=c(-.2, .2)) +
  scale_x_continuous(breaks=seq(0, 4, 1)) +
  theme(panel.border = element_blank(), panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"))+
  theme(axis.line.x = element_line(color="black", size = .5),
        axis.line.y = element_line(color="black", size = .5))+
  geom_hline(yintercept=0, color="gray")+
  xlab("Time From Adoption (in Years)") + ylab("Estimated Coefficient with CI")+
  ggtitle("Republican Ideology")+theme(plot.title = element_text(hjust = 0.5))

plot(plot.ideo.gop)

###############SLES

#keep only relevant variables, doesn't matter but panelmatch throws an error
dat.sles<-dat[c("sles", "party_id", "treatment", "majority", "party_id", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro", "year", "distid1")]

#have to constantly be making these as integers, super annoying
dat.sles<-as.data.frame(dat.sles)
dat.sles$distid1<-as.integer(dat.sles$distid1)
dat.sles$year<-as.integer(dat.sles$year)

# Create PanelData object for SLES
panel_data_sles <- PanelData(
  panel.data = dat.sles,
  unit.id = "distid1",
  time.id = "year",
  treatment = "treatment",
  outcome = "sles"
)

results.ps.sles <- PanelMatch(
  panel.data = panel_data_sles,
  lag = 4,
  refinement.method = "ps.match",
  match.missing = TRUE,
  covs.formula = ~ I(lag(majority, 1:4))+ I(lag(party_id, 1:4)) + I(lag(veto, 1:4)) + 
    I(lag(mds1, 1:4)) + I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4)) +
    I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
    I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
  size.match = 10,
  qoi = "att",
  lead = 0:4,
  forbid.treatment.reversal = TRUE
)

mset.sles <- results.ps.sles$att
PE.results.sles <- PanelEstimate(sets = results.ps.sles, panel.data = panel_data_sles)
summary(PE.results.sles)
PE.results.sles<-summary(PE.results.sles) #save as dataframe
names(PE.results.sles)
coef<-PE.results.sles[ , 1] #get time coefficients
se<-PE.results.sles[ , 2] #get standard errors
lower<-PE.results.sles[ , 3] #lower confidence bound
upper<-PE.results.sles[ , 4] #upper confidence bound
time<-c(0,1,2,3,4)

PE.results.sles<-as.data.frame(cbind(coef,se,lower,upper,time))
PE.results.sles

plot.sles<-ggplot(PE.results.sles, aes(time,coef)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) +
  theme_classic() +
  scale_y_continuous(breaks=seq(-1, 1, .4)) +
  coord_cartesian(ylim=c(-1, 1)) +
  scale_x_continuous(breaks=seq(0, 4, 1)) +
  theme(panel.border = element_blank(), panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"))+
  theme(axis.line.x = element_line(color="black", size = .5),
        axis.line.y = element_line(color="black", size = .5))+
  geom_hline(yintercept=0, color="gray")+
  xlab("Time From Adoption (in Years)") + ylab("Estimated Coefficient with CI")+
  ggtitle("Legislative Effectiveness")+theme(plot.title = element_text(hjust = 0.5))

plot(plot.sles)

###############Party Loyalty

#keep only relevant variables, doesn't matter but panelmatch throws an error
dat.loyalty<-dat[c("party_loyalty", "party_id", "treatment", "majority", "party_id", "veto", "mds1", "mds2", "state_ideo", "ranney", "lngsp", "tl", "lnexp", "citizen_ideo", "lnpop", "total_intro", "year", "distid1")]

#have to constantly be making these as integers, super annoying
dat.loyalty<-as.data.frame(dat.loyalty)
dat.loyalty$distid1<-as.integer(dat.loyalty$distid1)
dat.loyalty$year<-as.integer(dat.loyalty$year)

# Create PanelData object for Party Loyalty
panel_data_loyalty <- PanelData(
  panel.data = dat.loyalty,
  unit.id = "distid1",
  time.id = "year",
  treatment = "treatment",
  outcome = "party_loyalty"
)

results.psmatch.loyalty <- PanelMatch(
  panel.data = panel_data_loyalty,
  lag = 4,
  refinement.method = "ps.match",
  match.missing = TRUE,
  covs.formula = ~ I(lag(majority, 1:4))+ I(lag(party_id, 1:4)) + I(lag(veto, 1:4)) + 
    I(lag(mds1, 1:4)) + I(lag(mds2, 1:4)) + I(lag(state_ideo, 1:4)) +
    I(lag(ranney, 1:4)) + I(lag(lngsp, 1:4)) + I(lag(tl, 1:4)) + I(lag(lnexp, 1:4)) + 
    I(lag(citizen_ideo, 1:4)) + I(lag(lnpop, 1:4)) + I(lag(total_intro, 1:4)),
  size.match = 10,
  qoi = "att",
  lead = 0:4,
  forbid.treatment.reversal = TRUE
)

mset.loyalty <- results.psmatch.loyalty$att
PE.results.loyalty <- PanelEstimate(sets = results.psmatch.loyalty, panel.data = panel_data_loyalty)
summary(PE.results.loyalty)
PE.results.loyalty<-summary(PE.results.loyalty) #save as dataframe
names(PE.results.loyalty)
coef<-PE.results.loyalty[ , 1] #get time coefficients
se<-PE.results.loyalty[ , 2] #get standard errors
lower<-PE.results.loyalty[ , 3] #lower confidence bound
upper<-PE.results.loyalty[ , 4] #upper confidence bound
time<-c(0,1,2,3,4)

PE.results.loyalty<-as.data.frame(cbind(coef,se,lower,upper,time))
PE.results.loyalty

plot.loyalty<-ggplot(PE.results.loyalty, aes(time,coef)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) +
  theme_classic() +
  scale_y_continuous(breaks=seq(-.2, .2, .04)) +
  coord_cartesian(ylim=c(-.2, .2)) +
  scale_x_continuous(breaks=seq(0, 4, 1)) +
  theme(panel.border = element_blank(), panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.background = element_blank(),
        axis.line = element_line(colour = "black"))+
  theme(axis.line.x = element_line(color="black", size = .5),
        axis.line.y = element_line(color="black", size = .5))+
  geom_hline(yintercept=0, color="gray")+
  xlab("Time From Adoption (in Years)") + ylab("Estimated Coefficient with CI")+
  ggtitle("Party Loyalty")+theme(plot.title = element_text(hjust = 0.5))

plot(plot.loyalty)

all.plots<-grid.arrange(plot.ideo.dem, plot.ideo.gop, plot.sles, plot.loyalty, ncol=2)

plot(all.plots)

pdf("treat_effects_leg_level_v1.pdf", width = 14, height = 8.5)
plot(all.plots)
dev.off()


#Figure I5 - Updated for newer version
par(mfrow=c(2,3))
plot(mset.ideo.dem, xlim = c(0, 1000), ylim=c(0,250), main="Dem. Party Ideology")
plot(mset.ideo.gop, xlim = c(0, 1000),ylim=c(0,250), main="Rep. Party Ideology")
plot(mset.sles, xlim = c(0,1000),ylim=c(0,250),  main="Legislative Effectiveness")
plot(mset.loyalty, xlim = c(0, 1000),ylim=c(0,250), main="Party Loyalty")